library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(lattice)
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(sp)
library(rgdal)
## rgdal: version: 1.5-18, (SVN revision 1082)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: D:/Program Files/R-4.0.3/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: D:/Program Files/R-4.0.3/library/rgdal/proj
## Linking to sp version:1.4-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
library(tmap)
library(ggplot2)
library(gridExtra)
library(gstat)
library(OpenStreetMap)
library(spacetime)
library(knitr)
library(tmap)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## √ tibble 3.0.4 √ dplyr 1.0.2
## √ tidyr 1.1.2 √ stringr 1.4.0
## √ readr 1.4.0 √ forcats 0.5.0
## √ purrr 0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::combine() masks gridExtra::combine()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(spdep)
library(sf)
library(OpenStreetMap)
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: rpart
## Registered S3 method overwritten by 'spatstat':
## method from
## print.boxx cli
##
## spatstat 1.64-1 (nickname: 'Help you I can, yes!')
## For an introduction to spatstat, type 'beginner'
##
## Note: spatstat version 1.64-1 is out of date by more than a year; we strongly recommend upgrading to the latest version.
##
## Attaching package: 'spatstat'
## The following object is masked from 'package:gstat':
##
## idw
## The following object is masked from 'package:lattice':
##
## panel.histogram
library(here)
## here() starts at D:/CASA/STDM/practical 1
library(sp)
library(rgeos)
## rgeos version: 0.5-5, (SVN revision 640)
## GEOS runtime version: 3.8.0-CAPI-1.13.1
## Linking to sp version: 1.4-4
## Polygon checking: TRUE
library(maptools)
library(GISTools)
## Loading required package: RColorBrewer
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:spatstat':
##
## area
## The following object is masked from 'package:dplyr':
##
## select
library(tmap)
library(sf)
library(geojson)
##
## Attaching package: 'geojson'
## The following object is masked from 'package:graphics':
##
## polygon
library(geojsonio)
##
## Attaching package: 'geojsonio'
## The following object is masked from 'package:base':
##
## pretty
library(tmaptools)
dt <- st_read("data/factorscore.geojson")
## Reading layer `factorscore' from data source `D:\CASA\STDM\practical 1\Data\factorscore.geojson' using driver `GeoJSON'
## Simple feature collection with 779 features and 5 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 523850.3 ymin: 177849.8 xmax: 531172.2 ymax: 183893.7
## geographic CRS: WGS 84
dt <- dt %>%
st_set_crs(., 27700) #adjust the proj
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
##dt <- st_transform(dt, 4326)
# the spatial matrix
W_cont_el <- poly2nb(dt, queen=T)
#Queen's Case
W_cont_el_mat <- nb2listw(W_cont_el, style="W", zero.policy=TRUE)
#calculate Getis-Ord Gi* scores
lg1 <- localG(dt$region_1, listw=W_cont_el_mat, zero.policy=T)
dt$lg1 <- lg1[]
dt <- as(dt, 'Spatial')
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in CRS
## definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in CRS definition
tmaplondon <- dt %>%
st_bbox(.) %>%
tmaptools::read_osm(., type = "osm", zoom = NULL) ##london osm map
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
breaks1<-c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)
GIColours<- rev(brewer.pal(8, "RdBu"))
tmap_mode("plot")
## tmap mode set to plotting
#now plot on an interactive map
tm <- qtm(tmaplondon) + tm_shape(dt) +
tm_polygons("lg1",
style="fixed",
breaks=breaks1,
palette=GIColours,
midpoint=NA,
alpha = 0.7,
title="Getis-Ord Gi* of factor0") + tm_layout(legend.position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7)+
tm_compass(position = c("left", "bottom"),type = "arrow") +
tm_scale_bar(position = c("left", "bottom"))
#tmap_save(tm, filename = "map1.png")
tm
###############################
lg2 <- localG(dt$region_2, listw=W_cont_el_mat, zero.policy=T)
dt$lg2 <- lg2[]
tm1 <- qtm(tmaplondon) + tm_shape(dt) +
tm_polygons("lg2",
style="fixed",
breaks=breaks1,
palette=GIColours,
midpoint=NA,
alpha = 0.7,
title="Getis-Ord Gi* of factor1") + tm_layout(legend.position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7)+
tm_compass(position = c("left", "bottom"),type = "arrow") +
tm_scale_bar(position = c("left", "bottom"))
tm1
#tmap_save(tm1, filename = "map2.png")
#########################
lg3 <- localG(dt$region_3, listw=W_cont_el_mat, zero.policy=T)
dt$lg3 <- lg3[]
tm2 <- qtm(tmaplondon) + tm_shape(dt) +
tm_polygons("lg3",
style="fixed",
breaks=breaks1,
palette=GIColours,
midpoint=NA,
alpha = 0.7,
title="Getis-Ord Gi* of factor2") + tm_layout(legend.position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7)+
tm_compass(position = c("left", "bottom"),type = "arrow") +
tm_scale_bar(position = c("left", "bottom"))
tm2
#tmap_save(tm2, filename = "map3.png")
############################
dt$Kmeans_Cluster = as.character(dt$Kmeans_Cluster)
tmc <- qtm(tmaplondon)+tm_shape(dt) +
tm_polygons("Kmeans_Cluster", alpha = 0.85, palette = 'Pastel1', title="Result of Cluster", border.alpha = 0.5
)+ tm_layout(legend.position = c("right", "top"), legend.bg.color = "white", legend.bg.alpha = 0.7)+
tm_compass(position = c("left", "bottom"),type = "arrow") +
tm_scale_bar(position = c("left", "bottom"))
#tmap_save(tmc, filename = "mapc.png")
tmc
